# Title     : TODO
# Objective : TODO
# Created by: Administrator
# Created on: 2019/6/18

library(gplots)
library(Hmisc)
library(optparse)
library(ropls)
library(magrittr)
library(ComplexHeatmap)
library(tools)
library(ppcor)
library(tidyverse)

option_list <- list(
  make_option("--i", default = "AllMet_Raw.txt", type = "character", help = "metabolite data file"),
  make_option("--g", default = "SampleInfo.csv", type = "character", help = "sample group file"),
  make_option("--ef", default = "Extra_AllMet_Raw.txt", type = "character", help = "metabolite data file"),
  make_option("--confounder", default = "confounder_data.txt", type = "character", help = "confounder data file"),
  make_option("--config", default = "calculate_config.json", type = "character", help = "config file")
)
opt <- parse_args(OptionParser(option_list = option_list))

args <- commandArgs(trailingOnly = F)
scriptPath <- dirname(sub("--file=", "", args[grep("--file", args)]))
source(str_c(scriptPath, "/base.R"))

createWhenNoExist <- function(f) {
  !dir.exists(f) && dir.create(f)
}

corTest = function(da1, da2, method = "spearman") {
  name1 <- colnames(da1)
  name2 <- colnames(da2)
  tr_da1 <- as.matrix(da1)
  tr_da2 <- as.matrix(da2)
  num1 <- ncol(da1)
  num2 <- ncol(da2)
  pvalue <- vector()
  value <- vector()
  id1 <- vector()
  id2 <- vector()
  rec <- 1;
  for (i in 1:num1) {
    for (j in 1:num2) {
      a1 <- tr_da1[, i]
      a2 <- tr_da2[, j]
      id1[rec] <- name1[i]
      id2[rec] <- name2[j]
      a1 <- as.numeric(a1)
      a2 <- as.numeric(a2)
      corr <- cor.test(a1, a2, method = method)
      esti <- as.matrix(corr$estimate)[1]
      value[rec] <- esti
      pvalue[rec] <- corr$p.value
      rec <- rec + 1
    }
  }
  list(node1 = id1, node2 = id2, cor = value, p = pvalue)
}

configJson <- fromJSON(opt$config)
corP <- configJson$p %>%
  as.numeric()
corFdr <- configJson$fdr %>%
  as.numeric()
cor <- configJson$coe %>%
  as.numeric()
corMethod <- configJson$corMethod
fdrMethod <- configJson$fdrMethod
nodeSize <- configJson$nodeSize %>%
  as.numeric()
hasConfounder <- configJson$hasConfounder %>%
  as.logical()


orignalData <- read_tsv(opt$i) %>%
  as.data.frame()

allData <- if (!hasConfounder) {
  data <- orignalData %>%
    gather("SampleID", "Value", -Metabolite) %>%
    spread(Metabolite, "Value")
  parent <- "./"
  createWhenNoExist(parent)
  extraData <- read_tsv(opt$ef) %>%
    rename(Metabolite = 1) %>%
    gather("SampleID", "Value", -Metabolite) %>%
    spread(Metabolite, "Value")
  dataIds <- data$SampleID
  extraIds <- extraData$SampleID
  interIds <- intersect(dataIds, extraIds)
  extraData <- extraData %>%
    filter(SampleID %in% interIds) %>%
    arrange(factor(SampleID, levels = interIds)) %>%
    column_to_rownames("SampleID")
  inData <- data %>%
    filter(SampleID %in% interIds) %>%
    arrange(factor(SampleID, levels = interIds)) %>%
    column_to_rownames("SampleID")

  if (nrow(inData) < 2) {
    next
  }
  listRs <- corTest(inData, extraData, method = corMethod)
  data.frame(Node1 = listRs$node1, Node2 = listRs$node2, r = listRs$cor, P = listRs$p,
             stringsAsFactors = F)
}else {
  data <- orignalData %>%
    gather("ID", "Value", -Metabolite) %>%
    spread(Metabolite, "Value")
  extraData <- read_tsv(opt$ef) %>%
    rename(Metabolite = 1) %>%
    gather("ID", "Value", -Metabolite) %>%
    spread(Metabolite, "Value")
  confounderData <- read_tsv(opt$confounder) %>%
    rename(Metabolite = 1) %>%
    gather("ID", "Value", -Metabolite) %>%
    spread(Metabolite, "Value")
  dataIds <- data$ID
  extraIds <- extraData$ID
  confounderIds <- confounderData$ID
  interIds <- Reduce(intersect, list(dataIds, extraIds, confounderIds))
  confounderData <- confounderData %>%
    filter(ID %in% interIds) %>%
    arrange(factor(ID, levels = interIds))
  inData <- data %>%
    filter(ID %in% interIds) %>%
    arrange(factor(ID, levels = interIds))
  extraData <- extraData %>%
    filter(ID %in% interIds) %>%
    arrange(factor(ID, levels = interIds))
  if (nrow(inData) < 2) {
    next
  }
  pcorData <- zd_pcor(left = inData, right = extraData, confounding = confounderData, method = corMethod)
 tibble(Node1 = pcorData$pair_1, Node2 = pcorData$pair_2, r = pcorData$r, P = pcorData$p)
}

allData <- allData %>%
  mutate_at(vars(c("r")), function(x) {
    ifelse(is.na(x), 0, x)
  }) %>%
  mutate_at(vars(c("P")), function(x) {
    ifelse(is.na(x), 1, x)
  }) %>%
  mutate(FDR = p.adjust(P, method = fdrMethod))

corData <- allData %>%
  select(c("Node1", "Node2", "r")) %>%
  spread(Node1, "r") %>%
  mutate_at(vars(-c("Node2")), function(x) {
    ifelse(is.na(x), 0, x)
  }) %>%
  rename(` ` = Node2)
write.csv(corData, str_c("r_Matrix.csv"), quote = T, row.names = F)

pData <- allData %>%
  select(c("Node1", "Node2", "P")) %>%
  spread(Node1, "P") %>%
  rename(` ` = Node2)
write.csv(pData, str_c("P_Matrix.csv"), quote = T, row.names = F)

fdrData <- allData %>%
  select(c("Node1", "Node2", "FDR")) %>%
  spread(Node1, "FDR") %>%
  rename(` ` = Node2)
write.csv(fdrData, str_c("Corrected_P_Matrix.csv"), quote = T, row.names = F)

corFdr
edgeData <- allData %>%
  filter(P < corP & FDR < corFdr & abs(r) > cor) %>%
  filter(Node1 != Node2) %>%
  mutate(distName = {
    Node1 %>%
      map2_chr(Node2, function(x, y) {
        vec <- c(x, y) %>%
          sort()
        str_c(vec, collapse = ";")
      })
  }) %>%
  distinct(distName, .keep_all = T) %>%
  select(-c("distName")) %>%
  arrange(desc(abs(r)))
write_csv(edgeData, str_c("Network_Edges_for_Cytoscape.csv"))

nodes <- unique(c(edgeData$Node1, edgeData$Node2))
infoData <- tibble(Node = nodes) %>%
  mutate(Size = {
    Node %>%
      map_int(function(x) {
        edgeData %>%
          filter(Node1 == x | Node2 == x) %>%
          nrow()
      })
  }) %>%
  arrange(desc(Size)) %>%
  filter(Size >= nodeSize)
write_csv(infoData, str_c("Network_Nodes_for_Cytoscape.csv"))


